home *** CD-ROM | disk | FTP | other *** search
/ CD Exchange / CD Exchange - Volume 1.iso / education / airfoil / airfoil.c < prev    next >
C/C++ Source or Header  |  1989-09-24  |  13KB  |  489 lines
  1. /********************************************************************
  2. *                                                                   *
  3. * Joukowski Airfoil Generator with Streamline and Pressure          *
  4. * Distribution Algorithms                                           *
  5. *                                                                   *
  6. * Written by:   Russell Leighton      15 March 1987                 *
  7. *               Lancaster, CA                                       *
  8. * Modified by:  David Foster          19 June 1988                  *
  9. *               Rochester, MI                                       *
  10. *               To include effect of induced circulation            *
  11. *               by a first order approximation.                     *
  12. ********************************************************************/
  13.  
  14. #include "airfoil.h"
  15.  
  16. main()
  17. {
  18.    float rs,theta,h,vel,eta,S;
  19.    int psi;
  20.    ULONG MessageClass;
  21.  
  22.    open_things();
  23.    do_about();
  24.  
  25.    /****************/
  26.    /* Loop forever */
  27.    /****************/
  28.  
  29.    for(;;)
  30.    {
  31.       while (Continue)
  32.       {
  33.          /****************************************************/
  34.          /* Wait, initially and after each plot, for user to */
  35.          /* bring up the double menu requester               */
  36.          /****************************************************/
  37.  
  38.          Wait(1L<<w->UserPort->mp_SigBit);
  39.          if((message = (struct IntuiMessage *)GetMsg(w->UserPort)) != NULL)
  40.          {
  41.             MessageClass = message->Class;
  42.             ReplyMsg(message);
  43.             if (MessageClass == REQVERIFY)
  44.             {
  45.                do_request();
  46.                break;
  47.             }
  48.          } /* end if */
  49.       } /* end while */
  50.  
  51.       do_init();
  52.  
  53.       if(mode)
  54.       {
  55.  
  56.          /********************/
  57.          /* Plot Streamlines */
  58.          /********************/
  59.  
  60.          FILL = TRUE;
  61.          for (psi = 12;psi > 0;--psi)
  62.          {
  63.             do_mess();
  64.             if(!Continue) break;
  65.  
  66.             PENUP;
  67.             SetAPen(rp,(long)(psi+1));
  68.             SetBPen(rp,(long)(psi+1));
  69.  
  70.             for (theta = 0.015;theta < TWO_PI;theta += PI/100)
  71.             {
  72.                vel = psi/(4.*velocity*r*sin(theta));
  73.                S = (1.+sin(alpha)/sin(theta));
  74.                vel = abs(S/(2.*vel));
  75.                eta = sqrt(vel*vel+1.0)-vel;
  76.                rs = r*(1.+eta)/(1.-eta);
  77.                transform(rs,theta);
  78.             } /* end for */
  79.             AreaEnd(rp);
  80.          } /* end for */
  81.  
  82.          /*******************************/
  83.          /* Plot Stagnation Streamlines */
  84.          /*******************************/
  85.  
  86.          do_mess();
  87.          if(Continue)
  88.          {
  89.             FILL = FALSE;
  90.             PENUP;
  91.             SetAPen(rp,1L);
  92.             h = (r-4.0)/40.0;
  93.             theta = -alpha;
  94.  
  95.             for (rs = 4;rs >= r;rs += h)
  96.             {
  97.                transform(rs,theta);
  98.             }
  99.  
  100.             PENUP;
  101.             theta = PI + alpha;
  102.  
  103.             for (rs = 4;rs > r;rs += h)
  104.             {
  105.                transform(rs,theta);
  106.             }
  107.          } /* end if */
  108.       } /* end if */
  109.  
  110.       else
  111.       {
  112.  
  113.          /******************************/
  114.          /* Plot Pressure Distribution */
  115.          /******************************/
  116.  
  117.          do_mess();
  118.          if(Continue)
  119.          {
  120.             FILL = TRUE;
  121.             PENUP;
  122.             SetAPen(rp,2L);
  123.             SetBPen(rp,2L);
  124.  
  125.             for (theta = 0.0;theta <= TWO_PI;theta += PI/100)
  126.             {
  127.                rs = r+(sin(theta)+sin(alpha))*(sin(theta)+sin(alpha));
  128.                transform(rs,theta);
  129.             } /* end for */
  130.             AreaEnd(rp);
  131.          } /* end if */
  132.       } /* end else */
  133.  
  134.       /****************/
  135.       /* Plot Airfoil */
  136.       /****************/
  137.  
  138.       do_mess();
  139.       if(Continue)
  140.       {
  141.          FILL = TRUE;
  142.          PENUP;
  143.          rs = r;
  144.          SetAPen(rp,1L);
  145.          SetBPen(rp,1L);
  146.  
  147.          for (theta = 0.0;theta <= TWO_PI;theta += PI/100)
  148.          {
  149.             transform(rs,theta);
  150.          }
  151.          AreaEnd(rp);
  152.       } /* end if */
  153.    } /* end forever */
  154. } /* end main */
  155.  
  156. do_init()
  157. {
  158.    float a0;
  159.  
  160.    SetAPen(rp,0L);
  161.    RectFill(rp,(long)(XMIN+1),(long)(YMIN+1),(long)(XMAX-1),(long)(YMAX-1));
  162.    SetOPen(rp,1L);
  163.  
  164.    /***********************************************************/
  165.    /* Calculate circle constants (circle center and radius)   */
  166.    /* from airfoil constants through a reverse transformation */
  167.    /***********************************************************/
  168.  
  169.    alpha = (float)angle*PI/180;
  170.    c = (float)camber/25;
  171.    t = (float)thickness/25;
  172.    a = 0;
  173.    b = c/2;
  174.  
  175.    do
  176.    {
  177.       a0 = a;
  178.       a = t*(2*a0+1)/4/sqrt(b*b+2*a0+1);
  179.       b = c*(1+2*a)/(2+2*a);
  180.    }
  181.    while(abs(a-a0) > TOL);
  182.    
  183.    r = sqrt(b*b+(a+1)*(a+1));
  184.    Continue = TRUE;
  185. }
  186.  
  187. do_mess()
  188. {
  189.    ULONG MessageClass;
  190.  
  191.    /******************************************************************/
  192.    /* Check for double menu requester. This can be brought up at any */
  193.    /* time but can not be displayed unless the message is replied    */
  194.    /* to, therefore this routine must be called periodically.        */
  195.    /******************************************************************/
  196.  
  197.    if((message = (struct IntuiMessage *)GetMsg(w->UserPort)) != NULL)
  198.    {
  199.       MessageClass = message->Class;
  200.       ReplyMsg(message);
  201.       if(MessageClass == REQVERIFY) do_request();
  202.    }
  203.    ixo = XCEN;
  204.    iyo = YCEN;
  205. }
  206.  
  207. transform(rs,theta)
  208.  
  209. float rs,theta;
  210. {
  211.    float x,y,z,u,v;
  212.    int PLOT;
  213.    long ix,iy,cx,cy;
  214.  
  215.    /********************************************************************/
  216.    /* This is the Joukowski transformation routine. This is also where */
  217.    /* the plotting is done. Plotting is usually done using the area    */
  218.    /* fill routines (AreaDraw & AreaMove). If FILL is true then the    */
  219.    /* points are used to build up the area shape to be filled. See the */
  220.    /* Rom Kernal manual containing the graphics primatives for more    */
  221.    /* details. If FILL is false then normal plotting is done.          */
  222.    /********************************************************************/
  223.  
  224.    x = rs*cos(theta-alpha)+a;
  225.    y = rs*sin(theta-alpha)+b;
  226.    z = 1/(x*x+y*y);
  227.    u = x*(1+z)*cos(alpha)-y*(1-z)*sin(alpha);
  228.    v = y*(1-z)*cos(alpha)+x*(1+z)*sin(alpha);
  229.  
  230.    ix = (long)(SFAC*u+XCEN);
  231.    iy = (long)(YCEN-SFAC*v);
  232.  
  233.    if(FILL)
  234.    {
  235.       PLOT = FALSE;
  236.       cx = ix;
  237.       cy = iy;
  238.  
  239.       /*************************************************************/
  240.       /* This subroutine also clips the display area, hence the    */
  241.       /* following statements. Contrary to the what the Rom Kernal */
  242.       /* manual states, all clipping must be done by the program   */
  243.       /* if the area fill routines are used otherwise a nasty      */
  244.       /* crash will result if plotting takes place out of bounds.  */
  245.       /* Only the x values are clipped accurately since, for this  */
  246.       /* routine only the x values at the boundary need to be      */
  247.       /* accurate. The PEN parameter indicates the pen status for  */
  248.       /* moves and draws as set by either PENUP or PENDOWN.        */
  249.       /*************************************************************/
  250.  
  251.       if((ix <= XMIN) && (ixo > XMIN))
  252.       {
  253.          cy += (float)(iyo-iy)*(XMIN-ix)/(ixo-ix);
  254.          cx = XMIN;
  255.       }
  256.       else if((ixo <= XMIN) && (ix > XMIN))
  257.       {
  258.          iyo += (float)(iy-iyo)*(XMIN-ixo)/(ix-ixo);
  259.          ixo = XMIN;
  260.          AreaDraw(rp,ixo,iyo);
  261.          PLOT = TRUE;
  262.       }
  263.       else if((ix >= XMAX) && (ixo < XMAX))
  264.       {
  265.          cy += (float)(iyo-iy)*(XMAX-ix)/(ixo-ix);
  266.          cx = XMAX;
  267.       }
  268.       else if((ixo >= XMAX) && (ix < XMAX))
  269.       {
  270.          iyo += (float)(iy-iyo)*(XMAX-ixo)/(ix-ixo);
  271.          ixo = XMAX;
  272.          AreaDraw(rp,ixo,iyo);
  273.          PLOT = TRUE;
  274.       }
  275.  
  276.       if((ixo > XMIN) && (ixo < XMAX)) PLOT = TRUE;
  277.       if(cy < YMIN) cy = YMIN;
  278.       if(cy > YMAX) cy = YMAX;
  279.  
  280.       if(PLOT)
  281.       {
  282.          if(PEN) AreaDraw(rp,cx,cy);
  283.          else { AreaMove(rp,cx,cy); PENDOWN; }
  284.       }
  285.       ixo = ix;
  286.       iyo = iy;
  287.    }
  288.    else
  289.    {
  290.       if(PEN) Draw(rp,ix,iy);
  291.       else { Move(rp,ix,iy); PENDOWN; }
  292.    }
  293. }
  294.  
  295. do_request()
  296. {